#### SET UP ####
library(tidyverse) # For data cleaning
library(haven) # To read .sav files
library(psych) # For Cronbach's alpha
library(ggcorrplot) # For correlation matrix
library(survey) # For weighted stats
library(gtsummary) # For descriptive tables
library(performance) # To test if a model is good
library(stargazer) # For regression tables
library(ggeffects) # For interaction plots
library(Hmisc) # For data descriptives
library(broom) # For regression plots
library(patchwork) # For combining plots
library(sjPlot) # For interaction plots
library(openxlsx) # To save a table
library(lavaan) # For CFA
library(miceadds) # For clustered OLS
library(estimatr) # For clustered interaction plots
library(fabricatr) # For clustered interaction plots

# Read in citizen data
cit <- read_sav("2024-12-02cit_clean.sav") %>% 
  mutate(across(where(is.labelled), ~as_factor(.x)))

# Read in politician data
pol <- read_sav("2024-12-02pol_clean.sav") %>% 
  mutate(across(where(is.labelled), ~as_factor(.x)))


# Create dichotomous variable of political experience (> 0 -> experienced)
cit <- cit %>%
  mutate(pol_exp_dich = case_when(
    is.na(pol_exp_3) ~ NA_character_,
    pol_exp_3 > 0 ~ "Yes",
    TRUE ~ "No")) %>%
  mutate(pol_exp_dich = factor(pol_exp_dich, levels = c("No", "Yes")))

pol <- pol %>%
  mutate(pol_exp_dich = case_when(
    is.na(pol_exp) ~ NA_character_, 
    pol_exp > 0 ~ "Yes",
    TRUE ~ "No")) %>%
  mutate(pol_exp_dich = factor(pol_exp_dich, levels = c("No", "Yes"))) 

## Weight
## Citizen weight
cit %>% 
  summarise(
    min = min(rake_weights),
    max = max(rake_weights),
    mean = mean(rake_weights),
    median = median(rake_weights),
    sd = sd(rake_weights),
    IQR = IQR(rake_weights))

cit %>%
  ggplot(aes(x = rake_weights)) +
  geom_histogram(bins = 100) +
  labs(title = "", x = "Weights", y = "Count") +
  theme_bw()

ggsave("./Figures and tables/cit_wt_hist.png", dpi = 300,
       width = 8, height = 4)

## Politician weight
pol %>% 
  summarise(
    min = min(rake_weights),
    max = max(rake_weights),
    mean = mean(rake_weights),
    median = median(rake_weights),
    sd = sd(rake_weights),
    IQR = IQR(rake_weights))

pol %>%
  ggplot(aes(x = rake_weights)) +
  geom_histogram(bins = 100) +
  labs(title = "", x = "Weights", y = "Count") +
  theme_bw()

ggsave("./Figures and tables/pol_wt_hist.png", dpi = 300,
       width = 8, height = 4)


#### DEMOGRAPHICS ####
## Table AX

# Citizen - Unweighted
tbl <- cit %>% 
  tbl_summary(
    include = c(gender, agegrp4, edu2, domicil),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      gender = "Gender",
      agegrp4 = "Age group",
      edu2 = "University or university of applied sciences degree",
      domicil = "Type of community"),
    missing = "no") %>% 
  modify_header(label = "**Variable**")

tbl

tbl %>%
  as_hux_xlsx("./Figures and tables/demo_cit_unweighted.xlsx")

# Citizen - Weighted
tbl <- cit %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    include = c(gender, agegrp4, edu2, domicil),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      gender = "Gender",
      agegrp4 = "Age group",
      edu2 = "University or university of applied sciences degree",
      domicil = "Type of community"),
    missing = "no" ) %>% 
  modify_header(label = "**Variable**")

tbl

tbl %>%
  as_hux_xlsx("./Figures and tables/demo_cit_weighted.xlsx")

# Politicians - Unweighted
tbl <- pol %>% 
  tbl_summary(
    include = c(gender, agegrp4, edu2, domicil),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      gender = "Gender",
      agegrp4 = "Age group",
      edu2 = "University or university of applied sciences degree",
      domicil = "Type of community"),
    missing = "no") %>% 
  modify_header(label = "**Variable**")

tbl %>%
  as_hux_xlsx("./Figures and tables/demo_pol_unweighted.xlsx")

# Citizen - Weighted
tbl <- pol %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    include = c(gender, agegrp4, edu2, domicil),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      gender = "Gender",
      agegrp4 = "Age group",
      edu2 = "University or university of applied sciences degree",
      domicil = "Type of community"),
    missing = "no" ) %>% 
  modify_header(label = "**Variable**")

tbl %>%
  as_hux_xlsx("./Figures and tables/demo_pol_weighted.xlsx")




## Index
# EFA and Cronbach's alpha for the index for citizens
index <- cit %>% 
  select(multiculture:eu_influence) %>% 
  na.omit()

fa.parallel(index, fa = "fa")
KMO(cor(index))
cortest.bartlett(cor(index), n = 2689)
eigen(cor(index, use = "pairwise.complete.obs"))$values
fa(index, nfactors = 2, rotate = "oblimin", fm = "ml")
alpha(index, check.keys = TRUE)

# EFA and Cronbach's alpha for the index for politicians
index <- pol %>% 
  select(multiculture:eu_influence) %>% 
  na.omit()

fa.parallel(index, fa = "fa")
KMO(cor(index))
cortest.bartlett(cor(index), n = 2689)
eigen(cor(index, use = "pairwise.complete.obs"))$values
fa(index, nfactors = 2, rotate = "oblimin", fm = "ml")
alpha(index, check.keys = TRUE)

## Correlation matrix
# Citizens
corr_cit <- cit %>% 
  select(left_right, lib_con, sat_dem:illiberalism, rake_weights, -sat_dem_cat) %>% 
  na.omit()

corr_cit_wt <- cov.wt(corr_cit[, -11], wt = corr_cit$rake_weights, cor = TRUE)
cit_matrix <- corr_cit_wt$cor

ggcorrplot(cit_matrix,
           title = "",
           type = "lower",
           ggtheme = ggplot2::theme_bw,
           color = c("#440154FF", "#21908CFF", "#FDE725FF"),
           lab = TRUE,
           digits = 1)

ggsave("./Figures and tables/matrix_cit.png", dpi = 300,
       width = 8, height = 6)

# Politicians
corr_pol <- pol %>% 
  select(left_right, lib_con, sat_dem:illiberalism, rake_weights, -sat_dem_cat) %>% 
  na.omit()

corr_pol_wt <- cov.wt(corr_pol[, -11], wt = corr_pol$rake_weights, cor = TRUE)
cit_matrix <- corr_pol_wt$cor

ggcorrplot(cit_matrix,
           title = "",
           type = "lower",
           ggtheme = ggplot2::theme_bw,
           color = c("#440154FF", "#21908CFF", "#FDE725FF"),
           lab = TRUE,
           digits = 1)

ggsave("./Figures and tables/matrix_pol.png", dpi = 300,
       width = 8, height = 6)

#### DESCRIPTIVES ####

## Histograms and density plots of dependent variable
# Citizens

mean_illiberalism <- weighted.mean(cit$illiberalism, cit$rake_weights, na.rm = TRUE)
var_illiberalism <- wtd.var(cit$illiberalism, cit$rake_weights)
sd_illiberalism <- sqrt(var_illiberalism)

ggplot(cit, aes(x = illiberalism, weight = rake_weights)) +
  geom_histogram(aes(y = ..density..), bins = 50, fill = "lightblue", color = "black", alpha = 0.7) +
  geom_density(aes(y = ..density..), color = "black", size = 1, adjust = 1.5) +
  labs(title = "",
       x = "Illiberalism index",
       y = "Density") +
  annotate("text", x = 0.9, y = 2.5, 
           label = paste("Mean:", round(mean_illiberalism, 2), "\nSD:", round(sd_illiberalism, 2)), 
           hjust = 0, vjust = -5, size = 3, color = "black") +
  theme_bw()

ggsave("./Figures and tables/cit_index.png", dpi = 300,
       width = 8, height = 6)

# Politicians
mean_illiberalism <- weighted.mean(pol$illiberalism, pol$rake_weights, na.rm = TRUE)
var_illiberalism <- wtd.var(pol$illiberalism, pol$rake_weights)
sd_illiberalism <- sqrt(var_illiberalism)

ggplot(pol, aes(x = illiberalism, weight = rake_weights)) +
  geom_histogram(aes(y = ..density..), bins = 50, fill = "lightblue", color = "black", alpha = 0.7) +
  geom_density(aes(y = ..density..), color = "black", size = 1, adjust = 1.5) +
  labs(title = "",
       x = "Illiberalism index",
       y = "Density") +
  annotate("text", x = 0.9, y = 2.5, 
           label = paste("Mean:", round(mean_illiberalism, 2), "\nSD:", round(sd_illiberalism, 2)), 
           hjust = 0, vjust = -0.9, size = 3, color = "black") +
  theme_bw()

ggsave("./Figures and tables/pol_index.png", dpi = 300,
       width = 8, height = 6)


### Table 2
# Citizen - Weighted
tbl <- cit %>% 
  filter(!is.na(party_2)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_2,
    include = c(illiberalism, multiculture, women_rights, gay_rights, immigration, christian_values, eu_influence),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index",
      multiculture = "Multiculture",
      women_rights = "Women rights",
      gay_rights = "Gay rights",
      immigration = "Immigration",
      christian_values = "Christian values",
      eu_influence = "Globalism"),
    missing = "no" )  %>%
  add_overall(last = TRUE) %>% 
  add_p() %>% 
  modify_header(label = "**Variable**")

tbl %>%
  as_hux_xlsx("./Figures and tables/cit_index.xlsx")

# Politicians - Weighted
tbl <- pol %>% 
  filter(!is.na(party_2)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_2,
    include = c(illiberalism, multiculture, women_rights, gay_rights, immigration, christian_values, eu_influence),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index",
      multiculture = "Multiculture",
      women_rights = "Women rights",
      gay_rights = "Gay rights",
      immigration = "Immigration",
      christian_values = "Christian values",
      eu_influence = "Globalism"),
    missing = "no" )  %>% 
  add_overall(last = TRUE) %>% 
  add_p() %>% 
  modify_header(label = "**Variable**")

tbl %>%
  as_hux_xlsx("./Figures and tables/pol_index.xlsx")

### Table 3
cit %>% 
  filter(!is.na(party_2)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_2,
    include = c(illiberalism, multiculture, women_rights, gay_rights, immigration, christian_values, eu_influence),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index",
      multiculture = "Multiculture",
      women_rights = "Women rights",
      gay_rights = "Gay rights",
      immigration = "Immigration",
      christian_values = "Christian values",
      eu_influence = "Globalism"),
    missing = "no" )  %>%
  add_overall(last = TRUE) %>% 
  add_p() %>% 
  modify_header(label = "**Variable**")



pol %>% 
  filter(!is.na(party_2)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_2,
    include = c(illiberalism, multiculture, women_rights, gay_rights, immigration, christian_values, eu_influence),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index",
      multiculture = "Multiculture",
      women_rights = "Women rights",
      gay_rights = "Gay rights",
      immigration = "Immigration",
      christian_values = "Christian values",
      eu_influence = "Globalism"),
    missing = "no" )  %>%
  add_overall(last = TRUE) %>% 
  add_p() %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians")


#### Plot of table 1
variables <- c("illiberalism", "multiculture", "women_rights", "gay_rights", "immigration", "christian_values", "eu_influence")


pol_long <- pol %>%
  select(all_of(variables), rake_weights, party_2) %>%
  pivot_longer(cols = all_of(variables), names_to = "Variable", values_to = "Value") %>%
  mutate(Dataset = "pol")


cit_long <- cit %>%
  select(all_of(variables), rake_weights, party_2) %>%
  pivot_longer(cols = all_of(variables), names_to = "Variable", values_to = "Value") %>%
  mutate(Dataset = "cit")

# Combine the datasets
combined_data <- bind_rows(pol_long, cit_long)

combined_data <- combined_data %>% 
  filter(!is.na(Value) & !is.na(rake_weights) & !is.na(party_2)) %>% 
  mutate(
    Variable = factor(
      Variable,
      levels = c(
        "illiberalism", "multiculture", "women_rights", "gay_rights", 
        "immigration", "christian_values", "eu_influence"
      ),
      labels = c(
        "Illiberalism index", "Multiculturalism", "Gender equality", "Sexual minority rights",
        "Immigration", "Christian Values", "European integration"
      )
    ),
    Dataset = factor(
      Dataset,
      levels = c("cit", "pol"),
      labels = c("Voters", "Candidates")
    )
  )

mean_data <- combined_data %>%
  group_by(Variable, Dataset, party_2) %>%
  dplyr::summarize(
    Mean = weighted.mean(Value, rake_weights, na.rm = TRUE),
    .groups = "drop"
  )



mean_data <- mean_data %>%
  mutate(
    y_position = ifelse(Dataset == "Voters", 1, 1)
  )

## Free y-axis version
# ggplot(combined_data, aes(x = Value, weight = rake_weights, fill = Dataset)) +
#   geom_density(alpha = 0.5) +
#   facet_grid(Variable ~ party_2, scales = "free_y") +
#   geom_vline(
#     data = mean_data,
#     aes(xintercept = Mean, color = Dataset),
#     linetype = "dotted", size = 1,
#     inherit.aes = FALSE
#   ) +
#   geom_text(
#     data = mean_data,
#     aes(
#       x = Mean,
#       y = ifelse(Dataset == "Voters", 1, 1.5),
#       label = round(Mean, 3)
#     ),
#     hjust = -0.3, 
#     size = 3, 
#     fontface = "bold",
#     inherit.aes = FALSE
#   ) +
#   labs(
#     title = "",
#     x = "Illiberalism",
#     y = "Density",
#     fill = "",
#     color = ""
#   ) +
#   theme_bw() + 
#   theme(legend.position = "bottom")
# 
# 
# ggsave("./Figures and tables/attitudes_plot.png", dpi = 300,
#        width = 9, height = 10)



ggplot(combined_data, aes(x = Value, weight = rake_weights, fill = Dataset)) +
  geom_density(alpha = 0.5) +
  facet_grid(Variable ~ party_2) +
  geom_vline(
    data = mean_data,
    aes(xintercept = Mean, color = Dataset),
    linetype = "dotted", size = 1,
    inherit.aes = FALSE
  ) +
  geom_text(
    data = mean_data,
    aes(
      x = Mean,
      y = ifelse(Dataset == "Voters", 3, 4),
      label = round(Mean, 3)
    ),
    hjust = -0.3, 
    size = 3,
    fontface = "bold",
    inherit.aes = FALSE
  ) +
  labs(
    title = "",
    x = "Illiberalism",
    y = "Density",
    fill = "",
    color = ""
  ) +
  theme_bw() + 
  theme(legend.position = "bottom")


ggsave("./Figures and tables/attitudes_plot_not_free.png", dpi = 300,
       width = 9, height = 10)


## Table 4


# Tables for citizens
cit %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no" )  %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Citizens")

cit %>% 
  filter(!is.na(pol_exp_dich)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = pol_exp_dich,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>% 
  add_p() %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Citizens - Political experience")

cit %>% 
  filter(!is.na(sat_dem_cat)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = sat_dem_cat,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>% 
  add_p() %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Citizens - Satisfaction with democracy")

cit %>% 
  filter(!is.na(party_2)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_2,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>%
  add_p() %>%  
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Citizens - Party type")

cit %>% 
  filter(!is.na(left_right_cat)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = left_right_cat,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>%
  add_p() %>%  
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Citizens - Left-right self-identification")

cit %>% 
  filter(!is.na(lib_con_cat)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = lib_con_cat,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>%
  add_p() %>%  
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Citizens - Liberal-conservative self-identification")




# Tables for politicians
pol %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no" )  %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians")

pol %>% 
  filter(!is.na(pol_exp_dich)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = pol_exp_dich,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>% 
  add_p() %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians - Political experience")

pol %>% 
  filter(!is.na(sat_dem_cat)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = sat_dem_cat,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>% 
  add_p() %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians - Satisfaction with democracy")

pol %>% 
  filter(!is.na(party_2)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_2,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>%
  add_p() %>%  
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians - Party type")

pol %>% 
  filter(!is.na(left_right_cat)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = left_right_cat,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>%
  add_p() %>%  
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians - Left-right self-identification")

pol %>% 
  filter(!is.na(lib_con_cat)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = lib_con_cat,
    include = c(illiberalism),
    digits = list(all_continuous() ~ c(3, 3)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index"),
    missing = "no")  %>%
  add_p() %>%  
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians - Liberal-conservative self-identification")


### Table AX - Party wise
cit %>% 
  filter(!is.na(party_10)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_10,
    include = c(illiberalism, multiculture, women_rights, gay_rights, immigration, christian_values, eu_influence),
    digits = list(all_continuous() ~ c(2, 2)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index",
      multiculture = "Multiculture",
      women_rights = "Women rights",
      gay_rights = "Gay rights",
      immigration = "Immigration",
      christian_values = "Christian values",
      eu_influence = "Globalism"),
    missing = "no" )  %>%
  add_overall(last = TRUE) %>% 
  add_p() %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Citizens")

pol %>% 
  filter(!is.na(party_10)) %>% 
  svydesign(ids = ~1, data = ., weights = ~rake_weights) %>% 
  tbl_svysummary(
    by = party_10,
    include = c(illiberalism, multiculture, women_rights, gay_rights, immigration, christian_values, eu_influence),
    digits = list(all_continuous() ~ c(2, 2)),
    statistic = list(all_categorical() ~ "{n} ({p}%)",
                     all_continuous() ~ "{mean} ({sd})"),
    label = list(
      illiberalism = "Illiberalism index",
      multiculture = "Multiculture",
      women_rights = "Women rights",
      gay_rights = "Gay rights",
      immigration = "Immigration",
      christian_values = "Christian values",
      eu_influence = "Globalism"),
    missing = "no" )  %>%
  add_overall(last = TRUE) %>% 
  add_p() %>% 
  modify_header(label = "**Variable**") %>% 
  as_gt() %>%
  gt::tab_header(title = "Politicians")



### Figure of table 2

## Illiberalism index

pol_ill <- pol %>%
  select(illiberalism, rake_weights) %>%
  mutate(Dataset = "pol")

cit_ill <- cit %>%
  select(illiberalism, rake_weights) %>%
  mutate(Dataset = "cit")


# Combine the datasets
comb <- bind_rows(pol_ill, cit_ill)

comb <- comb %>% 
  filter(!is.na(illiberalism) & !is.na(rake_weights)) %>% 
  mutate(
    Dataset = factor(
      Dataset,
      levels = c("cit", "pol"),
      labels = c("Voters", "Candidates")
    )
  )

mean_data <- comb %>%
  group_by(Dataset) %>%
  dplyr::summarize(
    Mean = weighted.mean(illiberalism, rake_weights, na.rm = TRUE),
    .groups = "drop"
  )



mean_data <- mean_data %>%
  mutate(
    y_position = ifelse(Dataset == "Voters", 1, 1)
  )


p1 <- ggplot(comb, aes(x = illiberalism, weight = rake_weights, fill = Dataset)) +
  geom_density(alpha = 0.5) +
  #facet_grid( ~ Dataset) +
  geom_vline(
    data = mean_data,
    aes(xintercept = Mean, color = Dataset),
    linetype = "dotted", size = 1,
    inherit.aes = FALSE
  ) +
  geom_text(
    data = mean_data,
    aes(
      x = Mean,
      y = ifelse(Dataset == "Voters", 1, 2),
      label = round(Mean, 3)
    ),
    hjust = -0.3, 
    size = 3,
    fontface = "bold",
    inherit.aes = FALSE
  ) +
  labs(
    title = "Illiberalism index",
    x = "Illiberalism",
    y = "Density",
    fill = "",
    color = ""
  ) +
  theme_bw() +
  scale_y_continuous(limits = c(0, 4.3)) +
  theme(legend.position = "none")



## Political experience

pol_exp <- pol %>%
  select(illiberalism, pol_exp_dich, rake_weights) %>%
  mutate(Dataset = "pol")

cit_exp <- cit %>%
  select(illiberalism, pol_exp_dich, rake_weights) %>%
  mutate(Dataset = "cit")


# Combine the datasets
comb <- bind_rows(pol_exp, cit_exp)

comb <- comb %>% 
  filter(!is.na(illiberalism) & !is.na(rake_weights) & !is.na(pol_exp_dich)) %>% 
  mutate(
    Dataset = factor(
      Dataset,
      levels = c("cit", "pol"),
      labels = c("Voters", "Candidates")
    ),
    pol_exp_dich = factor(
      pol_exp_dich,
      levels = c("No", "Yes")
    )
  )

mean_data <- comb %>%
  group_by(Dataset, pol_exp_dich) %>%
  dplyr::summarize(
    Mean = weighted.mean(illiberalism, rake_weights, na.rm = TRUE),
    .groups = "drop"
  )



mean_data <- mean_data %>%
  mutate(
    y_position = ifelse(Dataset == "Voters", 1, 2)
  )


p2 <- ggplot(comb, aes(x = illiberalism, weight = rake_weights, fill = Dataset)) +
  geom_density(alpha = 0.5) +
  facet_grid( ~ pol_exp_dich) +
  geom_vline(
    data = mean_data,
    aes(xintercept = Mean, color = Dataset),
    linetype = "dotted", size = 1,
    inherit.aes = FALSE
  ) +
  geom_text(
    data = mean_data,
    aes(
      x = Mean,
      y = ifelse(Dataset == "Voters", 1, 1.5),
      label = round(Mean, 3)
    ),
    hjust = -0.3, 
    size = 3,
    fontface = "bold",
    inherit.aes = FALSE
  ) +
  labs(
    title = "Political experience",
    x = "Illiberalism",
    y = "Density",
    fill = "",
    color = ""
  ) +
  theme_bw() +
  scale_y_continuous(limits = c(0, 4.3))  +
  theme(legend.position = "none")


## Satisfaction with democracy

pol_sat <- pol %>%
  select(illiberalism, sat_dem, sat_dem_cat, rake_weights) %>%
  mutate(Dataset = "pol")

cit_sat <- cit %>%
  select(illiberalism, sat_dem, sat_dem_cat, rake_weights) %>%
  mutate(Dataset = "cit")

comb <- bind_rows(pol_sat, cit_sat)

comb <- comb %>% 
  filter(!is.na(illiberalism) & !is.na(rake_weights) & !is.na(sat_dem_cat)) %>% 
  mutate(
    Dataset = factor(
      Dataset,
      levels = c("cit", "pol"),
      labels = c("Voters", "Candidates")
    )
  )

mean_data <- comb %>%
  group_by(Dataset, sat_dem_cat) %>%
  dplyr::summarize(
    Mean = weighted.mean(illiberalism, rake_weights, na.rm = TRUE),
    .groups = "drop"
  )



mean_data <- mean_data %>%
  mutate(
    y_position = ifelse(Dataset == "Voters", 1, 1)
  )

p3 <- ggplot(comb, aes(x = illiberalism, weight = rake_weights, fill = Dataset)) +
  geom_density(alpha = 0.5) +
  facet_grid( ~sat_dem_cat) +
  geom_vline(
    data = mean_data,
    aes(xintercept = Mean, color = Dataset),
    linetype = "dotted", size = 1,
    inherit.aes = FALSE
  ) +
  geom_text(
    data = mean_data,
    aes(
      x = Mean,
      y = ifelse(Dataset == "Voters", 1, 2),
      label = round(Mean, 3)
    ),
    hjust = -0.3, 
    size = 3,
    fontface = "bold",
    inherit.aes = FALSE
  ) +
  labs(
    title = "Satisfaction with democracy",
    x = "Illiberalism",
    y = "Density",
    fill = "",
    color = ""
  ) +
  theme_bw() +
  scale_y_continuous(limits = c(0, 4.3)) +
  theme(legend.position = "none")  


## Party type

pol_comp <- pol %>%
  select(illiberalism, party_2, rake_weights) %>%
  mutate(Dataset = "pol")

cit_comp <- cit %>%
  select(illiberalism, party_2, rake_weights) %>%
  mutate(Dataset = "cit")

comb <- bind_rows(pol_comp, cit_comp)

comb <- comb %>% 
  filter(!is.na(illiberalism) & !is.na(rake_weights) & !is.na(party_2)) %>% 
  mutate(
    Dataset = factor(
      Dataset,
      levels = c("cit", "pol"),
      labels = c("Voters", "Candidates")
    )
  )

mean_data <- comb %>%
  group_by(Dataset, party_2) %>%
  dplyr::summarize(
    Mean = weighted.mean(illiberalism, rake_weights, na.rm = TRUE),
    .groups = "drop"
  )



mean_data <- mean_data %>%
  mutate(
    y_position = ifelse(Dataset == "Voters", 1, 1)
  )

p4 <- ggplot(comb, aes(x = illiberalism, weight = rake_weights, fill = Dataset)) +
  geom_density(alpha = 0.5) +
  facet_grid(~party_2 ) +
  geom_vline(
    data = mean_data,
    aes(xintercept = Mean, color = Dataset),
    linetype = "dotted", size = 1,
    inherit.aes = FALSE
  ) +
  geom_text(
    data = mean_data,
    aes(
      x = Mean,
      y = ifelse(Dataset == "Voters", 1, 2),
      label = round(Mean, 3)
    ),
    hjust = -0.3, 
    size = 3,
    fontface = "bold",
    inherit.aes = FALSE
  ) +
  labs(
    title = "Party type",
    x = "Illiberalism",
    y = "Density",
    fill = "",
    color = ""
  ) +
  theme_bw() +
  scale_y_continuous(limits = c(0, 4.3)) +
  theme(legend.position = "none")  



## Left-right

pol_comp <- pol %>%
  select(illiberalism, left_right_cat, rake_weights) %>%
  mutate(Dataset = "pol")

cit_comp <- cit %>%
  select(illiberalism, left_right_cat, rake_weights) %>%
  mutate(Dataset = "cit")

comb <- bind_rows(pol_comp, cit_comp)

comb <- comb %>% 
  filter(!is.na(illiberalism) & !is.na(rake_weights) & !is.na(left_right_cat)) %>% 
  mutate(
    Dataset = factor(
      Dataset,
      levels = c("cit", "pol"),
      labels = c("Voters", "Candidates")
    )
  )

mean_data <- comb %>%
  group_by(Dataset, left_right_cat) %>%
  dplyr::summarize(
    Mean = weighted.mean(illiberalism, rake_weights, na.rm = TRUE),
    .groups = "drop"
  )



mean_data <- mean_data %>%
  mutate(
    y_position = ifelse(Dataset == "Voters", 1, 1)
  )

p5 <- ggplot(comb, aes(x = illiberalism, weight = rake_weights, fill = Dataset)) +
  geom_density(alpha = 0.5) +
  facet_grid(~left_right_cat ) +
  geom_vline(
    data = mean_data,
    aes(xintercept = Mean, color = Dataset),
    linetype = "dotted", size = 1,
    inherit.aes = FALSE
  ) +
  geom_text(
    data = mean_data,
    aes(
      x = Mean,
      y = ifelse(Dataset == "Voters", 1, 2),
      label = round(Mean, 3)
    ),
    hjust = -0.3, 
    size = 3,
    fontface = "bold",
    inherit.aes = FALSE
  ) +
  labs(
    title = "Left-right self-identification",
    x = "Illiberalism",
    y = "Density",
    fill = "",
    color = ""
  ) +
  theme_bw() +
  scale_y_continuous(limits = c(0, 4.3)) +
  theme(legend.position = "none")  

## Liberal-conservative

pol_comp <- pol %>%
  select(illiberalism, lib_con_cat, rake_weights) %>%
  mutate(Dataset = "pol")

cit_comp <- cit %>%
  select(illiberalism, lib_con_cat, rake_weights) %>%
  mutate(Dataset = "cit")

comb <- bind_rows(pol_comp, cit_comp)

comb <- comb %>% 
  filter(!is.na(illiberalism) & !is.na(rake_weights) & !is.na(lib_con_cat)) %>% 
  mutate(
    Dataset = factor(
      Dataset,
      levels = c("cit", "pol"),
      labels = c("Voters", "Candidates")
    )
  )

mean_data <- comb %>%
  group_by(Dataset, lib_con_cat) %>%
  dplyr::summarize(
    Mean = weighted.mean(illiberalism, rake_weights, na.rm = TRUE),
    .groups = "drop"
  )



mean_data <- mean_data %>%
  mutate(
    y_position = ifelse(Dataset == "Voters", 1, 1)
  )

p6 <- ggplot(comb, aes(x = illiberalism, weight = rake_weights, fill = Dataset)) +
  geom_density(alpha = 0.5) +
  facet_grid(~lib_con_cat ) +
  geom_vline(
    data = mean_data,
    aes(xintercept = Mean, color = Dataset),
    linetype = "dotted", size = 1,
    inherit.aes = FALSE
  ) +
  geom_text(
    data = mean_data,
    aes(
      x = Mean,
      y = ifelse(Dataset == "Voters", 1, 2),
      label = round(Mean, 3)
    ),
    hjust = -0.3, 
    size = 3,
    fontface = "bold",
    inherit.aes = FALSE
  ) +
  labs(
    title = "Liberal-conservative self-identification",
    x = "Illiberalism",
    y = "Density",
    fill = "",
    color = ""
  ) +
  theme_bw() +
  scale_y_continuous(limits = c(0, 4.3)) +
  theme(legend.position = "bottom")  


p1 / p2 / p3 / p4 / p5 / p6



ggsave("./Figures and tables/attitudes_groupwise.png", dpi = 300,
       width = 9, height = 10)







#### REGRESSIONS ####


## Citizens
cit_reg <- cit %>% 
  filter(!is.na(party_10))

glimpse(cit)
cit %>% 
  count(party_10)

# Political experience
cit1_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(cit1_cl)
nobs(cit1_cl$lm_res)

# Satisfaction with democracy
cit2_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ sat_dem,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(cit2_cl)
nobs(cit2_cl$lm_res)

# Full attitudinal
cit3_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(cit3_cl)
nobs(cit3_cl$lm_res)

# Full attitudinal + controls
cit4_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(cit4_cl)
nobs(cit4_cl$lm_res)


## Run the same analyses as regular OLS -> see model fit and extract using stargazer
cit1 <- lm(illiberalism ~ pol_exp_dich, data = cit_reg, weights = rake_weights)
cit2 <- lm(illiberalism ~ sat_dem, data = cit_reg, weights = rake_weights)
cit3 <- lm(illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con, data = cit_reg, weights = rake_weights)
cit4 <- lm(illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens, data = cit_reg, weights = rake_weights)

stargazer(cit1, cit2, cit3, cit4, 
          type = "html",
          star.cutoffs = c(0.05, 0.01, 0.001),
          title = "Predicting illiberalism (citizens)",
          single.row=TRUE,
          dep.var.labels = c("Illiberalism index"), 
          covariate.labels = c("Political experience: Yes", "Satisfaction with democracy", "Party type: Illiberal", "Left-right identification", "Liberal-conservative identification", "Gender: Woman", "Age", "Education: University degree", "Population density")
)


# Model 4 best when comparing model performance indices
tbl <- compare_performance(cit1, cit2, cit3, cit4)

write.xlsx(tbl, file = "./Figures and tables/cit_fit.xlsx", rowNames = T)

# Check assumptions for the best model
check_model(cit4)


## Politicians
pol_reg <- pol %>% 
  filter(!is.na(party_10))

glimpse(pol)
pol %>% 
  count(party_10)

# Political experience
pol1_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(pol1_cl)
nobs(pol1_cl$lm_res)

# Satisfaction with democracy
pol2_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ sat_dem,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(pol2_cl)
nobs(pol2_cl$lm_res)

# Full attitudinal
pol3_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(pol3_cl)
nobs(pol3_cl$lm_res)

# Full attitudinal + controls
pol4_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(pol4_cl)
nobs(pol4_cl$lm_res)


## Run the same analyses as regular OLS -> see model fit and extract using stargazer
pol1 <- lm(illiberalism ~ pol_exp_dich, data = pol_reg, weights = rake_weights)
pol2 <- lm(illiberalism ~ sat_dem, data = pol_reg, weights = rake_weights)
pol3 <- lm(illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con, data = pol_reg, weights = rake_weights)
pol4 <- lm(illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens, data = pol_reg, weights = rake_weights)


stargazer(pol1, pol2, pol3, pol4, 
          type = "html",
          star.cutoffs = c(0.05, 0.01, 0.001),
          title = "Predicting illiberalism (politicians)",
          single.row=TRUE,
          dep.var.labels = c("Illiberalism index"), 
          covariate.labels = c("Political experience: Yes", "Satisfaction with democracy", "Party type: Illiberal", "Left-right identification", "Liberal-conservative identification", "Gender: Woman", "Age", "Education: University degree", "Population density")
)



# Model 4 best when comparing model performance indices
tbl <- compare_performance(pol1, pol2, pol3, pol4)

write.xlsx(tbl, file = "./Figures and tables/pol_fit.xlsx", rowNames = T)

# Check assumptions for the best model
check_model(pol4)



### Regression plots
## Voters

# Extract clustered CIs
sum_cit4_cl <- summary(cit4_cl)
coef_df <- as.data.frame(sum_cit4_cl)
coef_df <- cbind(term = rownames(coef_df), coef_df)

# Calculate confidence intervals
cit_r <- coef_df %>% 
  mutate(
    estimate = Estimate,
    conf.low = Estimate - 1.96 * `Std. Error`,
    conf.high = Estimate + 1.96 * `Std. Error`
  ) %>% 
  select(term, estimate, conf.low, conf.high)

# Define reference categories
reference_cat <- data.frame(
  term = c(
    "Political experience: No",
    "Party: Liberal",
    "Gender: Male",
    "Education: Doesn't have a degree"),
  estimate = c(0, 0, 0, 0),
  conf.low = c(0, 0, 0, 0),
  conf.high = c(0, 0, 0, 0))

cit_r <- cit_r %>% 
  filter(term != "(Intercept)") %>% 
  mutate(term = case_when(
    term == "pol_exp_dichYes" ~ "Political experience: Yes",
    term == "sat_dem" ~ "Satisfaction with democracy",
    term == "party_2Illiberal" ~ "Party: Illiberal",
    term == "left_right" ~ "Left-right identification",
    term == "lib_con" ~ "Liberal-conservative identification",
    term == "genderWoman" ~ "Gender: Female",
    term == "age" ~ "Age",
    term == "edu2Degree" ~ "Education: Has a degree",
    term == "pop_dens" ~ "Population density")) %>% 
  select(term, estimate, conf.low, conf.high)

cit_r <- rbind(cit_r, reference_cat)

cit_r <- cit_r %>% 
  mutate(term = factor(term, levels = c(
    "Population density",
    "Education: Has a degree",
    "Education: Doesn't have a degree",
    "Age",
    "Gender: Female",
    "Gender: Male",
    "Liberal-conservative identification",
    "Left-right identification",
    "Party: Illiberal",
    "Party: Liberal",
    "Satisfaction with democracy",
    "Political experience: Yes",
    "Political experience: No"
  )),
  type = "Voters")

# cit_r_p <- ggplot(cit_r, aes(x = term, y = estimate)) +
#   geom_hline(yintercept = 0, linetype = "dotted", color = "red") +
#   geom_point(size = 3) +
#   geom_linerange(aes(ymin = conf.low, ymax = conf.high)) +
#   labs(
#     title = "Voters",
#     x = "",
#     y = ""
#   ) +
#   coord_flip() +
#   theme_bw() +
#   scale_y_continuous(limits = c(-0.25, 0.5), expand = c(0, 0))






## Candidates
# Extract clustered CIs
sum_pol4_cl <- summary(pol4_cl)
coef_df <- as.data.frame(sum_pol4_cl)
coef_df <- cbind(term = rownames(coef_df), coef_df)

# Calculate confidence intervals
pol_r <- coef_df %>% 
  mutate(
    estimate = Estimate,
    conf.low = Estimate - 1.96 * `Std. Error`,
    conf.high = Estimate + 1.96 * `Std. Error`
  ) %>% 
  select(term, estimate, conf.low, conf.high)


pol_r <- pol_r %>% 
  filter(term != "(Intercept)") %>% 
  mutate(term = case_when(
    term == "pol_exp_dichYes" ~ "Political experience: Yes",
    term == "sat_dem" ~ "Satisfaction with democracy",
    term == "party_2Illiberal" ~ "Party: Illiberal",
    term == "left_right" ~ "Left-right identification",
    term == "lib_con" ~ "Liberal-conservative identification",
    term == "genderWoman" ~ "Gender: Female",
    term == "age" ~ "Age",
    term == "edu2Degree" ~ "Education: Has a degree",
    term == "pop_dens" ~ "Population density")) %>% 
  select(term, estimate, conf.low, conf.high)

pol_r <- rbind(pol_r, reference_cat)

pol_r <- pol_r %>% 
  mutate(term = factor(term, levels = c(
    "Population density",
    "Education: Has a degree",
    "Education: Doesn't have a degree",
    "Age",
    "Gender: Female",
    "Gender: Male",
    "Liberal-conservative identification",
    "Left-right identification",
    "Party: Illiberal",
    "Party: Liberal",
    "Satisfaction with democracy",
    "Political experience: Yes",
    "Political experience: No"
  )),
  type = "Candidates")

# pol_r_p <- ggplot(pol_r, aes(x = term, y = estimate)) +
#   geom_hline(yintercept = 0, linetype = "dotted", color = "red") +
#   geom_point(size = 3) +
#   geom_linerange(aes(ymin = conf.low, ymax = conf.high)) +
#   labs(
#     title = "Candidates",
#     x = "",
#     y = ""
#   ) +
#   coord_flip() +
#   theme_bw() +
#   scale_y_continuous(limits = c(-0.25, 0.5), expand = c(0, 0))
# 
# cit_r_p / pol_r_p


## Combined
comb <- rbind(cit_r, pol_r)

comb <- as.data.frame(comb)

comb <- comb %>% 
  mutate(type = factor(type, levels = c("Voters", "Candidates")))

ggplot(comb, aes(x = term, y = estimate, color = type, shape = type)) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), 
                 position = position_dodge(width = 0.5), size = 1) +
  labs(
    title = "",
    x = "",
    y = "",
    shape = "",
    color = ""
  ) +
  coord_flip() +
  theme_bw() +
  #scale_y_continuous(limits = c(-0.3, 0.6), expand = c(0, 0)) +
  theme(legend.position = "bottom")  

ggsave("./Figures and tables/pred_illiberalism.png", dpi = 300, 
       width = 8, height = 6)


#### INTERACTIONS ####

# Citizens - Interactions

citint1 <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich * sat_dem,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(citint1)
nobs(citint1$lm_res)

citint2 <- miceadds::lm.cluster(
  formula = illiberalism ~ party_2 * sat_dem,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(citint2)
nobs(citint2$lm_res)

citint3 <- miceadds::lm.cluster(
  formula = illiberalism ~ party_2 * pol_exp_dich,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(citint3)
nobs(citint3$lm_res)



# Politicians - Interactions

polint1 <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich * sat_dem,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(polint1)
nobs(polint1$lm_res)

polint2 <- miceadds::lm.cluster(
  formula = illiberalism ~ party_2 * sat_dem,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(polint2)
nobs(polint2$lm_res)

polint3 <- miceadds::lm.cluster(
  formula = illiberalism ~ party_2 * pol_exp_dich,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(polint3)
nobs(polint3$lm_res)



# Satdem - polexp

mod <- lm_robust(illiberalism ~  sat_dem * pol_exp_dich,
                 data    = cit_reg,
                 cluster = cit_reg$party_10,
                 weights = cit_reg$rake_weights,
                 se = "stata"
)

int1a <- plot_model(mod, type = "int",
           terms = c("sat_dem", "pol_exp_dich"), 
           title = "Voters") + 
  theme_bw() +
  labs(
    x = "Satisfaction with Democracy",
    y = "Illiberalism index",
    color = "Political experience",  
    fill = "Political experience"    
  )


mod <- lm_robust(illiberalism ~  sat_dem * pol_exp_dich,
                 data    = pol_reg,
                 cluster = pol_reg$party_10,
                 weights = pol_reg$rake_weights,
                 se = "stata"
)

int1b <- plot_model(mod, type = "int",
                    terms = c("sat_dem", "pol_exp_dich"), 
                    title = "Candidates") + 
  theme_bw() +
  labs(
    x = "Satisfaction with Democracy",
    y = "Illiberalism index",
    color = "Political experience",  
    fill = "Political experience"    
  )


int1a / int1b

ggsave("./Figures and tables/int_satdem_polexp.png", 
       dpi = 300, 
       width = 8, 
       height = 6)

# party2 - satdem

mod <- lm_robust(illiberalism ~  sat_dem * party_2,
                 data    = cit_reg,
                 cluster = cit_reg$party_10,
                 weights = cit_reg$rake_weights,
                 se = "stata"
)

int2a <- plot_model(mod, type = "int",
                    terms = c("sat_dem", "party_2"), 
                    title = "Voters") + 
  theme_bw() +
  labs(
    x = "Satisfaction with Democracy",
    y = "Illiberalism index",
    color = "Party type",  
    fill = "Party type"    
  )


mod <- lm_robust(illiberalism ~  sat_dem * party_2,
                 data    = pol_reg,
                 cluster = pol_reg$party_10,
                 weights = pol_reg$rake_weights,
                 se = "stata"
)

int2b <- plot_model(mod, type = "int",
                    terms = c("sat_dem", "party_2"), 
                    title = "Candidates") + 
  theme_bw() +
  labs(
    x = "Satisfaction with Democracy",
    y = "Illiberalism index",
    color = "Party type",  
    fill = "Party type"    
  )


int2a / int2b


ggsave("./Figures and tables/int_satdem_party.png", 
       dpi = 300, 
       width = 8, 
       height = 6)

# party2 - polexp


mod <- lm_robust(illiberalism ~  pol_exp_dich * party_2,
                 data    = cit_reg,
                 cluster = cit_reg$party_10,
                 weights = cit_reg$rake_weights,
                 se = "stata"
)

int3a <- plot_model(mod, type = "int",
                    terms = c("pol_exp_dich", "party_2"), 
                    title = "Voters") + 
  theme_bw() +
  labs(
    x = "Political experience",
    y = "Illiberalism index",
    color = "Party type",  
    fill = "Party type"    
  )


mod <- lm_robust(illiberalism ~  pol_exp_dich * party_2,
                 data    = pol_reg,
                 cluster = pol_reg$party_10,
                 weights = pol_reg$rake_weights,
                 se = "stata"
)

int3b <- plot_model(mod, type = "int",
                    terms = c("pol_exp_dich", "party_2"), 
                    title = "Candidates") + 
  theme_bw() +
  labs(
    x = "Political experience",
    y = "Illiberalism index",
    color = "Party type",  
    fill = "Party type"    
  )


int3a / int3b



ggsave("./Figures and tables/int_polexp_party.png", 
       dpi = 300, 
       width = 8, 
       height = 6)


#### ROBUSTNESS CHECKS ####

## Political experience as numeric

### Political experience citizens

exp1 <- lm(illiberalism ~ pol_exp_1, data = cit, weights = rake_weights)
exp2 <- lm(illiberalism ~ pol_exp_2, data = cit, weights = rake_weights)
exp3 <- lm(illiberalism ~ pol_exp_3, data = cit, weights = rake_weights)
exp4 <- lm(illiberalism ~ pol_exp_4, data = cit, weights = rake_weights)
exp5 <- lm(illiberalism ~ pol_exp_5, data = cit, weights = rake_weights)
exp6 <- lm(illiberalism ~ pol_exp_6, data = cit, weights = rake_weights)

tbl <- compare_performance(exp1, exp2, exp3, exp4, exp5, exp6)

write.xlsx(tbl, file = "./Figures and tables/cit_pol_exp_fit.xlsx", rowNames = T)



cit4 <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_3 + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens,
  data    = cit_reg,
  cluster = cit_reg$party_10,
  weights = cit_reg$rake_weights
)
summary(cit4)
nobs(cit4$lm_res)


pol4 <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_lvl + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens,
  data    = pol_reg,
  cluster = pol_reg$party_10,
  weights = pol_reg$rake_weights
)
summary(pol4)
nobs(pol4$lm_res)


### Regression plots


## Voters

# Extract clustered CIs
sum_cit4_cl <- summary(cit4)
coef_df <- as.data.frame(sum_cit4_cl)
coef_df <- cbind(term = rownames(coef_df), coef_df)

# Calculate confidence intervals
cit_r <- coef_df %>% 
  mutate(
    estimate = Estimate,
    conf.low = Estimate - 1.96 * `Std. Error`,
    conf.high = Estimate + 1.96 * `Std. Error`
  ) %>% 
  select(term, estimate, conf.low, conf.high)

# Define reference categories
reference_cat <- data.frame(
  term = c(
    "Party: Liberal",
    "Gender: Male",
    "Education: Doesn't have a degree"),
  estimate = c(0, 0, 0),
  conf.low = c(0, 0, 0),
  conf.high = c(0, 0, 0))

cit_r <- cit_r %>% 
  filter(term != "(Intercept)") %>% 
  mutate(term = case_when(
    term == "pol_exp_3" ~ "Political experience",
    term == "sat_dem" ~ "Satisfaction with democracy",
    term == "party_2Illiberal" ~ "Party: Illiberal",
    term == "left_right" ~ "Left-right identification",
    term == "lib_con" ~ "Liberal-conservative identification",
    term == "genderWoman" ~ "Gender: Female",
    term == "age" ~ "Age",
    term == "edu2Degree" ~ "Education: Has a degree",
    term == "pop_dens" ~ "Population density")) %>% 
  select(term, estimate, conf.low, conf.high)

cit_r <- rbind(cit_r, reference_cat)

cit_r <- cit_r %>% 
  mutate(term = factor(term, levels = c(
    "Population density",
    "Education: Has a degree",
    "Education: Doesn't have a degree",
    "Age",
    "Gender: Female",
    "Gender: Male",
    "Liberal-conservative identification",
    "Left-right identification",
    "Party: Illiberal",
    "Party: Liberal",
    "Satisfaction with democracy",
    "Political experience"
  )),
  type = "Voters")






## Candidates
# Extract clustered CIs
sum_pol4_cl <- summary(pol4)
coef_df <- as.data.frame(sum_pol4_cl)
coef_df <- cbind(term = rownames(coef_df), coef_df)

# Calculate confidence intervals
pol_r <- coef_df %>% 
  mutate(
    estimate = Estimate,
    conf.low = Estimate - 1.96 * `Std. Error`,
    conf.high = Estimate + 1.96 * `Std. Error`
  ) %>% 
  select(term, estimate, conf.low, conf.high)


pol_r <- pol_r %>% 
  filter(term != "(Intercept)") %>% 
  mutate(term = case_when(
    term == "pol_exp_lvl" ~ "Political experience",
    term == "sat_dem" ~ "Satisfaction with democracy",
    term == "party_2Illiberal" ~ "Party: Illiberal",
    term == "left_right" ~ "Left-right identification",
    term == "lib_con" ~ "Liberal-conservative identification",
    term == "genderWoman" ~ "Gender: Female",
    term == "age" ~ "Age",
    term == "edu2Degree" ~ "Education: Has a degree",
    term == "pop_dens" ~ "Population density")) %>% 
  select(term, estimate, conf.low, conf.high)

pol_r <- rbind(pol_r, reference_cat)

pol_r <- pol_r %>% 
  mutate(term = factor(term, levels = c(
    "Population density",
    "Education: Has a degree",
    "Education: Doesn't have a degree",
    "Age",
    "Gender: Female",
    "Gender: Male",
    "Liberal-conservative identification",
    "Left-right identification",
    "Party: Illiberal",
    "Party: Liberal",
    "Satisfaction with democracy",
    "Political experience"
  )),
  type = "Candidates")



## Combined
comb <- rbind(cit_r, pol_r)

comb <- as.data.frame(comb)

comb <- comb %>% 
  mutate(type = factor(type, levels = c("Voters", "Candidates")))

ggplot(comb, aes(x = term, y = estimate, color = type, shape = type)) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), 
                 position = position_dodge(width = 0.5), size = 1) +
  labs(
    title = "",
    x = "",
    y = "",
    shape = "",
    color = ""
  ) +
  coord_flip() +
  theme_bw() +
  #scale_y_continuous(limits = c(-0.3, 0.6), expand = c(0, 0)) +
  theme(legend.position = "bottom")  


ggsave("./Figures and tables/pol_exp_cont.png", dpi = 300, 
       width = 8, height = 6)


## Excluding "Other" parties
## Exclude those from parties not represented in the parliament

cit_oth <- cit_reg %>% 
  mutate(party_2 = case_when(
    party_10 == "Other" ~ NA_character_,
    TRUE ~ party_2
  )) %>% 
  mutate(party_2 = factor(party_2, levels = c("Liberal", "Illiberal")))

cit1_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich,
  data    = cit_oth,
  cluster = cit_oth$party_10,
  weights = cit_oth$rake_weights
)
summary(cit1_cl)
nobs(cit1_cl$lm_res)

cit2_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ sat_dem,
  data    = cit_oth,
  cluster = cit_oth$party_10,
  weights = cit_oth$rake_weights
)
summary(cit2_cl)
nobs(cit2_cl$lm_res)

cit3_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con,
  data    = cit_oth,
  cluster = cit_oth$party_10,
  weights = cit_oth$rake_weights
)
summary(cit3_cl)
nobs(cit3_cl$lm_res)

cit4_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens,
  data    = cit_oth,
  cluster = cit_oth$party_10,
  weights = cit_oth$rake_weights
)
summary(cit4_cl)
nobs(cit4_cl$lm_res)



# Politicians

pol_oth <- pol_reg %>% 
  mutate(party_2 = case_when(
    party_10 == "Other" ~ NA_character_,
    TRUE ~ party_2
  )) %>% 
  mutate(party_2 = factor(party_2, levels = c("Liberal", "Illiberal")))


pol1_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich,
  data    = pol_oth,
  cluster = pol_oth$party_10,
  weights = pol_oth$rake_weights
)
summary(pol1_cl)
nobs(pol1_cl$lm_res)

pol2_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ sat_dem,
  data    = pol_oth,
  cluster = pol_oth$party_10,
  weights = pol_oth$rake_weights
)
summary(pol2_cl)
nobs(pol2_cl$lm_res)

pol3_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con,
  data    = pol_oth,
  cluster = pol_oth$party_10,
  weights = pol_oth$rake_weights
)
summary(pol3_cl)
nobs(pol3_cl$lm_res)

pol4_cl <- miceadds::lm.cluster(
  formula = illiberalism ~ pol_exp_dich + sat_dem + party_2 + left_right + lib_con + gender + age + edu2 + pop_dens,
  data    = pol_oth,
  cluster = pol_oth$party_10,
  weights = pol_oth$rake_weights
)
summary(pol4_cl)
nobs(pol4_cl$lm_res)


#### BLACKBOX ####

## Set up

library(tidyverse) # For data cleaning
library(purrr) # For data cleaning
library(haven) # To read .sav files
library(basicspace) # For blackbox
library(asmcjr) # For model plotting, see here for installation: https://uniofessex.github.io/asmcjr/articles/installation.html
library(patchwork) # For combining plots

# https://uniofessex.github.io/asmcjr/articles/ch2.html

# Read in citizen data and create dichotomous political experience and sat dem variable
cit <- read_sav("2024-12-02cit_clean.sav") %>% 
  mutate(across(where(is.labelled), ~as_factor(.x)),
         type = "Voter") %>% 
  mutate(pol_exp_dich = case_when(
    is.na(pol_exp_3) ~ NA_character_,
    pol_exp_3 > 0 ~ "Yes",
    TRUE ~ "No")) %>%
  mutate(pol_exp_dich = factor(pol_exp_dich, levels = c("No", "Yes")),
         sat_dem_dich = if_else(sat_dem >= mean(sat_dem, na.rm = TRUE),
                                "High satisfaction", "Low satisfaction")
  )


cit %>% 
  count(sat_dem, sat_dem_dich)

glimpse(cit)

# Read in politician data and create dichotomous political experience variable
pol <- read_sav("2024-12-02pol_clean.sav") %>% 
  mutate(across(where(is.labelled), ~as_factor(.x)),
         type = "Candidate") %>%
  mutate(pol_exp_dich = case_when(
    is.na(pol_exp) ~ NA_character_, 
    pol_exp > 0 ~ "Yes",
    TRUE ~ "No")) %>%
  mutate(pol_exp_dich = factor(pol_exp_dich, levels = c("No", "Yes")),
         sat_dem_dich = if_else(sat_dem >= mean(sat_dem, na.rm = TRUE),
                                "High satisfaction", "Low satisfaction")
  ) 


# Combine the two datasets, only select relevant variables, do cleaning
comb <- bind_rows(cit, pol) %>% 
  select(multiculture:eu_influence, type, party_2, pol_exp_dich, sat_dem_dich) %>% 
  mutate(
    across(where(is.numeric), ~ as.integer(.x * 10)), # Multiply by 10 and make into integers
    across(where(is.numeric), ~ replace_na(., 999)), # Make NA into 999
    across(where(is.numeric), ~ as.numeric(as.character(.)))
  ) %>% 
  filter(!is.na(type), !is.na(party_2), !is.na(pol_exp_dich), !is.na(sat_dem_dich))

# Matrix with issues 
issues <- comb %>% 
  select(multiculture:eu_influence) %>% # Respondent type and individual index items
  as.matrix() 


# Run blackbox with two dimensions, minimum 3 responses 
two_dim <- blackbox(issues,
                    missing = 999,
                    dims = 2,
                    minscale = 3,
                    verbose = TRUE)
two_dim$fits




# Printing the R squared values of the issues by dimension, 
# using the two-dimensional estimation
two_dim$stimuli


# Type: voter or candidate
type <- comb %>% 
  select(type) %>% 
  mutate(type = factor(type))

# Party: liberal or illiberal
party <- comb %>% 
  select(party_2)

# Political experience: Yes or no
pol_exp <- comb %>% 
  select(pol_exp_dich)

# Satisfaction with democracy
sat_dem <- comb %>% 
  select(sat_dem_dich) %>% 
  mutate(sat_dem_dich = factor(sat_dem_dich))


# Plot two-dimensional blackbox by voters and candidates
p1 <- plot_blackbox(two_dim, 
                    dims=c(1,2),
                    groupVar = type$type,
                    main = "Respondent Type",
                    xlab = "First Dimension (58%)",
                    ylab ="Second Dimension (14%)") +
  labs(title = "Respondent Type") +
  theme(legend.position = "bottom", aspect.ratio=1,
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold")) +
  guides(shape = guide_legend(override.aes = list(size=4))) +
  scale_color_manual(values = c(
    "Candidate" = "#0072B2",
    "Voter" = "#D55E00"))

# Plot two-dimensional blackbox by party type
p2 <- plot_blackbox(two_dim, 
                    dims = c(1,2),
                    groupVar = party$party_2,
                    main = "Paty Type",
                    xlab = "First Dimension (58%)",
                    ylab = "Second Dimension (14%)") +
  labs(title = "Party Type") +
  theme(legend.position = "bottom", aspect.ratio=1,
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold")) +
  guides(shape = guide_legend(override.aes = list(size = 4))) +
  scale_color_manual(values = c(
    "Illiberal" = "#0072B2",
    "Liberal" = "#D55E00"))

# Plot two-dimensional blackbox by political experience
p3 <- plot_blackbox(two_dim, 
                    dims = c(1,2),
                    groupVar = pol_exp$pol_exp_dich,
                    xlab = "First Dimension (58%)",
                    ylab = "Second Dimension (14%)") +
  labs(title = "Political Experience") +
  theme(legend.position = "bottom", aspect.ratio=1,
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold")) +
  guides(shape = guide_legend(override.aes = list(size = 4))) +
  scale_color_manual(values = c(
    "No" = "#0072B2",
    "Yes" = "#D55E00"))

# Plot two-dimensional blackbox by dichotomous satisfaction with democracy
p4 <- plot_blackbox(two_dim, 
                    dims = c(1,2),
                    groupVar = sat_dem$sat_dem_dich,
                    xlab = "First Dimension (58%)",
                    ylab = "Second Dimension (14%)") +
  labs(title = "Satisfaction with Democracy") +
  theme(legend.position = "bottom", aspect.ratio=1,
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold")) +
  guides(shape = guide_legend(override.aes = list(size = 4))) +
  scale_color_manual(values = c(
    "High satisfaction" = "#0072B2",
    "Low satisfaction" = "#D55E00"))


(p1 + p2) / (p3 + p4)

ggsave("./Figures and tables/DIF/blackbox.png", 
       dpi = 300, 
       width = 8, 
       height = 8)


#### MULTI-GROUP CFA ####

glimpse(comb)

## CFA model
model <- 'illiberalism  =~ multiculture + immigration + women_rights + gay_rights + christian_values + eu_influence'

# Here tried to get the latent means of the groups, didn't really work 
# model <- 'illiberalism  =~ multiculture + immigration + women_rights + gay_rights + christian_values + eu_influence
#           illiberalism ~c(0,a)*1'

# Regular fit
fit <- cfa(model, data = comb, meanstructure = T)
summary(fit, standardized = T)

# CFI and TLI good, bad RMSEA, good srmr
fit_baseline <- fitMeasures(fit, c("cfi","tli","rmsea","srmr")) %>%
  as.data.frame() %>%
  mutate(model = "Baseline")
fit_baseline


## Want:
# CFI >= 0.9
# TLI >= 0.9
# RMSEA <= 0.08
# SRMR <= 0.08

## Multi-group CFA

## TYPE

## Configural invariance: same structure, different parameters
fit_configural <- cfa(model,
                      data = comb,
                      group = "type")
summary(fit_configural, fit.measures = TRUE)
fit_config_1 <- fitMeasures(fit_configural, c("cfi","tli","rmsea","srmr")) %>%
  as.data.frame() %>%
  mutate(model = "Configural (Type)")
fit_config_1

## Metric invariance: equal loadings
fit_metric <- cfa(model,
                  data = comb,
                  group = "type",
                  group.equal = "loadings")
summary(fit_metric, fit.measures = TRUE)

## Scalar invariance: equal loadings + intercepts
fit_scalar <- cfa(model,
                  data = comb,
                  group = "type",
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, fit.measures = TRUE)

tbl1 <- anova(fit_configural, fit_metric, fit_scalar)

tbl1
fitMeasures(fit_configural)

## Look at CFI, does it drop more than 0.01? Good or bad 
map_df(list(
    fit_configural, fit_metric, fit_scalar),
    function(x) fitMeasures(x, c("chisq", "df", "cfi", "rmsea", "srmr")))

## PARTY

## Configural invariance: same structure, different parameters
fit_configural <- cfa(model,
                      data = comb,
                      group = "party_2")
summary(fit_configural, fit.measures = TRUE)
fit_config_2 <- fitMeasures(fit_configural, c("cfi","tli","rmsea","srmr")) %>%
  as.data.frame() %>%
  mutate(model = "Configural (Party)")
fit_config_2

## Metric invariance: equal loadings
fit_metric <- cfa(model,
                  data = comb,
                  group = "party_2",
                  group.equal = "loadings")
summary(fit_metric, fit.measures = TRUE)

## Scalar invariance: equal loadings + intercepts
fit_scalar <- cfa(model,
                  data = comb,
                  group = "party_2",
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, fit.measures = TRUE)

tbl2 <- anova(fit_configural, fit_metric, fit_scalar)

tbl2


## EXPERIENCE

## Configural invariance: same structure, different parameters
fit_configural <- cfa(model,
                      data = comb,
                      group = "pol_exp_dich")
summary(fit_configural, fit.measures = TRUE)
fit_config_3 <- fitMeasures(fit_configural, c("cfi","tli","rmsea","srmr")) %>%
  as.data.frame() %>%
  mutate(model = "Configural (Experience)")
fit_config_3


## Metric invariance: equal loadings
fit_metric <- cfa(model,
                  data = comb,
                  group = "pol_exp_dich",
                  group.equal = "loadings")
summary(fit_metric, fit.measures = TRUE)

## Scalar invariance: equal loadings + intercepts
fit_scalar <- cfa(model,
                  data = comb,
                  group = "pol_exp_dich",
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, fit.measures = TRUE)

tbl3 <- anova(fit_configural, fit_metric, fit_scalar)

tbl3


## SATISFACTION

## Configural invariance: same structure, different parameters
fit_configural <- cfa(model,
                      data = comb,
                      group = "sat_dem_dich")
summary(fit_configural, fit.measures = TRUE)
fit_config_4 <- fitMeasures(fit_configural, c("cfi","tli","rmsea","srmr")) %>%
  as.data.frame() %>%
  mutate(model = "Configural (Satisfaction)")
fit_config_4

## Metric invariance: equal loadings
fit_metric <- cfa(model,
                  data = comb,
                  group = "sat_dem_dich",
                  group.equal = "loadings")
summary(fit_metric, fit.measures = TRUE)

## Scalar invariance: equal loadings + intercepts
fit_scalar <- cfa(model,
                  data = comb,
                  group = "sat_dem_dich",
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, fit.measures = TRUE)

tbl4 <- anova(fit_configural, fit_metric, fit_scalar)

tbl4


tbl <- bind_rows(tbl1, tbl2, tbl3, tbl4)

write.xlsx(tbl, file = "./Figures and tables/DIF/measurement_invariance.xlsx", rowNames = T)



tbl <- bind_rows(fit_baseline, fit_config_1, fit_config_2, fit_config_3, fit_config_4)
write.xlsx(tbl, file = "./Figures and tables/DIF/baseline_configural_fits.xlsx", rowNames = T)
